{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimation of gradient index n2 from measured focal range of GRIN-rod F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sympy as sym\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.rcParams.update({\n",
    "    'font.size': 16,\n",
    "    'legend.frameon': True\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Direct calculation of F from known n2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full formula\n",
    "\n",
    "The formula is initialy taken from the Yariv's 'Quantun electronics' p.113."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n0 = 1.73\n",
    "n2 = np.linspace(0.1, 100, 100)\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "ax.set_xlabel('Gradient constant $n_2 \\;\\; (1/m^2)$')\n",
    "ax.set_ylabel('Focus range, m')\n",
    "ax.set_yscale('log')\n",
    "ax.set_xscale('log')\n",
    "ax.grid(True, linestyle=':')\n",
    "for L in [0.001, 0.01, 0.1]:\n",
    "    F = 1/n0*np.sqrt(n0/n2)*np.cos(np.sqrt(n2/n0)*L)/np.sin(np.sqrt(n2/n0)*L)\n",
    "    ax.plot(n2, F, label=f'L={L}')\n",
    "ax.legend(title='Rod length, m')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Precision of approximation\n",
    "\n",
    "Approximation is done by means that `sin(x) = x` and `cos(x) = 1` when `x` is small."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n2 = np.linspace(0.1, 100, 10)\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "ax.set_xlabel('Gradient constant $n_2 \\;\\; (1/m^2)$')\n",
    "ax.set_ylabel('Approximation error, %')\n",
    "ax.grid(True, linestyle=':')\n",
    "for L in [0.01, 0.05, 0.1]:\n",
    "    F1 = 1/n0*np.sqrt(n0/n2)*np.cos(np.sqrt(n2/n0)*L)/np.sin(np.sqrt(n2/n0)*L)\n",
    "    F2 = 1/L/n2\n",
    "    ax.plot(n2, (F2-F1)/F1*100, label=f'L={L}')\n",
    "ax.legend(title='Rod length, m')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inverse calculation of n2 from measured F "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 0.1\n",
    "F = 1.45\n",
    "\n",
    "def equat_float(x):\n",
    "    return x*F*n0*math.tan(x) - L\n",
    "\n",
    "def equat_np(x):\n",
    "    return x*F*n0*np.tan(x) - L\n",
    "    \n",
    "def calc_F(n2):\n",
    "    return 1/n0*math.sqrt(n0/n2)/math.tan(math.sqrt(n2/n0)*L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculation from reduced formula\n",
    "\n",
    "x = sym.Symbol('x')\n",
    "x1 = [sym.nsolve(x*F*n0*sym.sin(x) - L*sym.cos(x), x0) for x0 in [0, 0.1, 0.5, 1]]\n",
    "print('Possible solutions for x:', x1)\n",
    "x_root = np.min([x for x in x1 if x > 0])\n",
    "n2_root = n0*(x_root/L)**2\n",
    "F1 = calc_F(n2_root)\n",
    "print('Solution: x_root =', x_root, 'n2 =', n2_root)\n",
    "print(f'Back calculated F={F1}, error={(F1-F)/F*100}%')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculation from full formula\n",
    "\n",
    "x = sym.Symbol('x')\n",
    "x1 = [sym.nsolve(sym.sqrt(n0/x)/n0*sym.cos(sym.sqrt(x/n0)*L) - F*sym.sin(sym.sqrt(x/n0)*L), x0) for x0 in [0.1, 0.5, 1]]\n",
    "print('Possible solutions for n2:', x1)\n",
    "n2 = np.min([x for x in x1 if x > 0])\n",
    "F1 = calc_F(n2)\n",
    "print('Solution: n2 =', n2)\n",
    "print(f'Diff with prev solution: diff_n2={abs(n2-n2_root)}')\n",
    "print(f'Back calculated F={F1}, error={(F1-F)/F*100}%')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate manually\n",
    "\n",
    "x = np.linspace(0, 0.3, 100)\n",
    "y = equat_np_neg(x)\n",
    "print('y_min =', np.min(y))\n",
    "print('y @ x_root', equat_float(x_root))\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "ax.grid(True, linestyle=':')\n",
    "ax.set_xlabel('$x = \\sqrt {\\\\frac {n_2} {n_0}} L $')\n",
    "ax.set_ylabel('$ x \\; F \\; n_0 \\; \\\\tan(x) - L $')\n",
    "ax.plot(x, y)\n",
    "ax.plot([np.min(x), np.max(x)], [0, 0])\n",
    "ax.plot([x_root, x_root], [np.min(y), np.max(y)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_next = 0\n",
    "x_prev = 0\n",
    "y_next = equat_float(x_next)\n",
    "iters1 = 0\n",
    "while y_next < 0:\n",
    "    iters1 += 1\n",
    "    x_prev = x_next\n",
    "    x_next += 0.1\n",
    "    y_next = equat_float(x_next)\n",
    "#print(f'x_prev={x_prev}, x_next={x_next}, y_next={y_next}')\n",
    "y_mid = y_next\n",
    "iters2 = 0\n",
    "while abs(y_mid) > 1e-8:\n",
    "    iters2 += 1\n",
    "    x_mid = (x_prev + x_next) / 2\n",
    "    y_mid = equat_float(x_mid)\n",
    "    #print(f'x_mid={x_mid}, y_mid={y_mid}, y_next={y_next}')\n",
    "    if y_next > 0 and y_mid < 0:\n",
    "        x_prev = x_mid\n",
    "    else:\n",
    "        x_next = x_mid\n",
    "        y_next = y_mid\n",
    "n2 = n0*(x_mid/L)**2\n",
    "F1 = calc_F(n2)\n",
    "print(f'Solution: x_mid={x_mid}, y_mid={y_mid}, n2={n2}, iters1={iters1}, iters2={iters2}, iters={iters1+iters2}')\n",
    "print(f'Back calculated F={F1}, error={abs(F1-F)/F*100}%')\n",
    "print(f'Diff with sym solution: diff_x={abs(x_mid-x_root)}, diff_n2={abs(n2-n2_root)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_r = np.linspace(-1, 1, 100)\n",
    "y_n = n0 - (n2*r**2)/2\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "ax.set_xlabel('Radial distance, m')\n",
    "ax.set_ylabel('IOR')\n",
    "ax.grid(True, linestyle=':')\n",
    "ax.plot(x_r, y_n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L_rod = L\n",
    "n2_this = n2\n",
    "x_n0 = np.linspace(1, 2.5, 100)\n",
    "y_F = 1/x_n0*np.sqrt(x_n0/n2)/np.tan(np.sqrt(n2/x_n0)*L)\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "ax.set_xlabel('$n_0$')\n",
    "ax.set_ylabel('Focal range, m')\n",
    "ax.grid(True, linestyle=':')\n",
    "ax.plot(x_n0, y_F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'n0={n0}, n2={n2}')\n",
    "x_L = np.linspace(0.01, 1, 100)\n",
    "y_F = 1/x_n0*np.sqrt(x_n0/n2)/np.tan(np.sqrt(n2/x_n0)*L)\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "ax.set_xlabel('Rod length, m')\n",
    "ax.set_ylabel('Focal range, m')\n",
    "ax.grid(True, linestyle=':')\n",
    "ax.plot(x_L, y_F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
